home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / FDMSTART / EX5 / CONSLAW.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  1.4 KB  |  52 lines
  1. #include <ConsLaw.h>
  2. #include <DrawFD.h>  // defines makeCurvPlot (FieldFD&,...)
  3. real start1 (const Ptv(real)& x, real t) { return 0.0; }
  4. real linearFlux (real u) { return u; }
  5.  
  6. void ConsLaw:: init ()
  7. {
  8.   f = linearFlux;
  9.  
  10.   grid.rebind (new GridLattice(1));  grid->scan (s_i);
  11.   u.rebind (new FieldFD (grid(), "u"));  
  12.   u_prev.rebind (new FieldFD (grid(), "u_prev"));
  13.   u_prev() = start1;
  14.   const real uL = 1;  u->valueIndex(0) = u_prev->valueIndex(0) = uL;
  15.   tip.scan (s_i);
  16.   file.open (CaseName);
  17.   time_points_for_plot.scan (s_i);
  18.  
  19.   // write input data for a check:
  20.   u->grid().print(s_o);  tip.print(s_o);
  21.   s_o << "time points for plot: "; time_points_for_plot.print (s_o);
  22. }
  23.  
  24. void ConsLaw:: solveProblem ()
  25. {
  26.   init ();
  27.  
  28.   // abbreviations:
  29. # define Un(i)  u->valueIndex(i)
  30. # define Un_1(i)  u_prev->valueIndex(i)
  31. # define Fn_1(i)  f(u_prev->valueIndex(i))
  32.  
  33.   real mu = tip.Delta()/u->grid().Delta(1);
  34.   const int n0 = u->grid().getBase(1);  // first grid point
  35.   const int n1 = u->grid().getMaxI(1);  // last  grid point
  36.   int i;
  37.   tip.initTimeLoop();
  38.   while (!tip.finished()) {
  39.     tip.increaseTime();
  40.     for (i = n0+1; i <= n1; i++)  {
  41.       Un(i) = Un_1(i) - mu * (Fn_1(i) - Fn_1(i-1));  }
  42.       if (time_points_for_plot.within (tip.getTime(), 0.4999*tip.Delta()))
  43.         DrawFD::makeCurvPlot (u(),file,"u(x) for t fixed","u",
  44.                   oform("t=%g",tip.getTime()));
  45.       u_prev() = u();
  46.     }
  47. # undef Un
  48. # undef Un_1
  49. # undef Fn_1
  50. }
  51.  
  52.